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摘 要 为 了 进一步 提高 辐射 传递 因子 的 倒 易 性 ， 本 文 提出 了 和 迭代 双向 统计 蒙特 卡 罗 方 法 . 该 方法 将 双向 统计 蒙特 
卡 罗 方 法 分 解 为 归 一 化 过 程 与 双向 统计 过 程 ， 摆 脱 了 传统 迭代 双向 统计 蒙特 卡 罗 方 法 对 等 权 抽样 假设 的 依赖 ， 进 一 
步 提高 了 辐射 传递 因子 的 倒 易 性 。 验 证 结 果 表明 ， 和 迭代 双向 统计 喧 特 卡 罗 方 法 获得 的 辐射 传递 因子 满足 完整 性 ， 其 
倒 易 性 误差 的 量 级 随 双向 统计 迁 代 次 数 的 增加 而 呈 线 性 下 降 趋 势 .将 该 方法 应 用 到 含 空 腔 结构 的 辐射 /导热 耦合 边界 
平衡 问题 中 , 结果 表明 本 文 所 建 方法 可 以 有 效 提高 热 辐射 传输 问题 的 求解 精度 , 并 且 求 解 精度 随 着 迭代 次 数 的 增加 而 提高 . 
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Iterative Bidirectional Statistics Monte Carlo Method in Thermal Radiation 
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Abstract An iterative bidirectional statistical Monte Carlo method is proposed to improve the 
reciprocity of radiative transfer factor, which decomposes the bidirectional statistical Monte Carlo 
method into a normalized process and a bidirectional statistical process. This method can get rid 
of the assumption of equal weight sampling and further improve the reciprocity of radiative transfer 
factor. The results of the validation of the new method show that the magnitude of the error of 
the reciprocity is decreased linearly with the increase of the number of iterations, and the integrity 
is certified. The method is applied into equilibrium radiation problem. The result shows that the 
method developed in this paper can improve the accuracy of the solution of radiation transfer, and 


the accuracy increase with the increase of the number of iterations. 
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0 引 言 


蒙特 卡 罗 方 法 是 一 种 概率 统计 方法 ， 既 能 用 于 
求解 随机 性 问题 ， 也 能 用 于 求解 确定 性 问题 .该 方 
法 最 早 由 Howell 和 Perlmutter 应 用 于 热 辐射 换 热 计 
算 143, 因 其 具有 物理 模型 简单 清晰 、 应 用 灵活 等 优 
点 , 获得 了 广泛 的 应 用 B- 引 , 但 是 , 由 于 该 方法 是 一 
种 统计 方法 ， 需 要 进行 大 量 的 抽样 模拟 以 提高 其 精 
度 . 一 些 学 者 结合 实际 问题 ， 提 出 了 一 些 可 以 改善 
其 计算 效率 的 处 理 方法 671， 其 中 就 包括 文献 [9] 
提出 的 双向 统计 蒙特 卡 罗 方 法 . 

双向 统计 蒙特 卡 罗 方 法 (BSMC) 采用 等 温 等 权 
抽样 假设 ,利用 热 辐 射 传输 时 光路 可 逆 原 理 ， 对 热 
辐射 传输 的 蒙特 卡 罗 模 拟 进行 双向 统计 ， 充 分 利用 
能 束 传播 的 路 径 信 息 ， 一定 程 度 上 提高 了 辐射 传递 
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因子 的 倒 易 性 四 . 但 该 方法 得 到 的 辐射 传递 因子 的 
倒 易 性 有 时 候 并 不 一 定 能 满足 实际 需求 . 为 了 更 进 
一 步 地 利用 能 束 传播 路 径 信息 ， 提 高 辐射 传递 因子 
的 倒 易 性 ， 本 文 提 出 了 迭代 双向 统计 蒙特 卡 罗 方 法 
(IBSMC)， 


1 友 代 双 加 统计 驼 特 卡 罗 方 法 


1.1 现 有 双向 统计 蒙特 卡 罗 方 法 
辐射 传递 因子 RDi 定义 为 在 一 个 红外 辐射 传 
输 系统 中 ,单元 i 的 本 身 辐射 能 量 经 系统 各 单元 一 
次 或 多 次 反射 和 散射 后 最 终 被 单元 7 吸收 的 份额 
以 由 MM 个 单元 构成 的 辐射 传输 系统 为 例 . 假设 
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系统 模拟 能 束 总 量 为 N， 单元 i 模拟 发 射 的 能 束 量 
为 Ni, 则 


i=M 
>》 N=N (1) 
| 
传统 蒙特 卡 罗 方 法 得 到 辐射 传递 因子 记 为 
RD9,, 则 
RD = 台 @) 
其 中 , Ni; 为 单元 i 发 射 的 能 束 最 终 被 单元 7 吸收 的 
能 束 量 ， 
双向 统计 蒙特 卡 罗 方 法 得 到 的 辐射 传递 因子 记 
为 RD;, 则 


x Ni + Nji 
RD = — i (3) 
Ni 十 >》， Nei 
k=1 
其 中 , NN; 满足 本 
OO (4) 
N; Ek; 
对 于 面 元 
Eb; 一 cji4; (5) 
对 于 体 元 


si 为 面 元 i 的 平均 发 射 率 ，A; 为 面 元 i 的 面积 ,ki 
为 体 元 i 的 平均 发 射 率 , Vi 为 体 元 i 的 面积 ， 

式 (4) 为 等 权 抽 样 假设 的 表述 ， 即 假设 每 束 光 
束 携带 相同 的 能 量 。 
1.2 迭代 双向 统计 蒙特 卡 罗 方 法 

将 双向 统计 蒙特 卡 罗 方 法 分 解 为 归 一 化 过 程 与 
双向 统计 过 程 。 归 一 化 过 程 即 传统 蒙特 卡 罗 方 法 过 
程 , 如 式 (2) 所 示 . 双向 统计 过 程 视 为 对 传统 蒙特 卡 
罗 方 法 得 到 的 辐射 传递 因子 的 后 处 理 过 程 , 如 式 (7) 


所 示 . 
Ei;RD?; + E;RD;, 


RDi; = 一 人 (7) 
Ei+ 》 ERDY, 
= 


由 式 (2) 和 式 (7) 可 以 看 出 , 归 一 化 过 程 和 双向 
统计 过 程 都 不 依赖 于 等 权 抽 样 假设 ， 即 将 式 (3) 分 
解 为 归 一 化 过 程 和 双向 统计 过 程 后 ， 使 双向 统计 蒙 
特 卡 罗 方 法 摆脱 了 等 权 抽 样 假设 ( 式 4) 的 限制 , 增 
大 了 其 适用 范围 。 

双向 统计 过 程 利用 了 光路 可 道 原理 ， 相 当 于 将 
辐射 换 热 系统 内 的 抽样 能 束 总 数目 增加 一 倍 ， 相 应 
地 ， 每 一 能 束 的 能 量 减 半 外， 而 蒙特 卡 罗 方 法 的 精 


度 随 着 抽样 能 束 的 增加 而 提高 。 因 此 可 以 考虑 利用 
迭代 的 方法 进行 多 次 双向 统计 过 程 ， 则 第 n 次 双向 
统计 得 到 的 辐射 传递 因子 表述 为 
多 一 工 多 一 工 
RD = 人 (8) 
Eit+ >》 ERD?! 
k=1 
其 中 ，RD% 为 第 ”次 双向 统计 得 到 的 辐射 传递 因 
子 , n 三 1,2,3,…。 式 (2) 和 式 (8) 联合 起 来 即 为 大 
代 双 向 统计 蒙特 卡 罗 方 法 , 当 n = 1 时 , 方法 退化 
为 文献 [9] 中 的 双向 统计 蒙特 卡 罗 方 法 . 
为 验证 迭代 式 (8) 得 到 的 辐射 传递 因子 RD8 
的 完整 性 ， 先 假设 RD% 满足 完整 性 ， 则 对 于 i = 
1,2,3,… ;, M;n 二 1,2,3,…, 有 


BRD? 3 EiRDY + EBRD 
2 k=M 
条 二 >》 ERDY! 
Bs1 
1 一 AM 


>》， (EiRDY 十 瓦 RD1 
一 


k=M (9) 
Eit+ 》, ERD?! 
k=1 
l=M 1 一 AI 
EE RD® +T》 ERD 1 
l=0 l=0 一 1 
k=M 
Ei+ >》 ERDr! 
k=1 
时 1 一 AI 
>》 RD =1 (10) 
l=0 


由 式 (9)、(10) 可 知 RDS 满足 完整 性 的 前 提 是 
RD2 满足 完整 性 . 已 知 传统 蒙特 卡 罗 方 法 得 到 加 
射 传递 因子 的 初始 值 RDS 满足 完整 性 , 由 递 推 关系 
知 ，RDY 满足 完整 性 . 
2 方法 验证 

选取 如 图 1 所 示 正方 形 空 腔 ， 空 腔 表 面 均 为 浊 
灰 体 表面 ,发射 率 为 0.5, 每 个 表面 划分 为 一 个 网 格 
单元 ， 每 个 单元 模拟 光束 量 为 1 x 10。 倒 易 性 误差 
定义 为 
a 2 (Ei;RDi; — E; RD;i) 


| 
BRD BRD 关 1000| OY) 


个 易 性 误差 和 完整 性 误差 随 迭 代 次 数 的 变化 分 


别 如 图 2 和 图 3 所 示 ， 本 文 提 出 的 IBSMC 方法 计 
算得 到 的 辐射 传递 因子 倒 易 性 误差 的 量 级 随 双 向 统 
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计 和 迭代 次 数 的 增加 呈 线 性 下 降 趋 势 ， 双 向 统计 的 迭 
代 过 程 中 辐射 传递 因子 完整 性 保持 非常 好 ,证 明了 
IBSMC 方法 的 可 靠 性 . 


图 1 正方 形 空 腔 
Fig. 1 Square cavity 


倒 易 性 误差 


0 2 4 6 8 
双向 统计 迭代 次 数 
图 2 倒 易 性 误差 
Fig. 2 Error of the reciprocity 
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利用 传统 的 蒙特 卡 罗 方 法 计算 耗 时 7.4 x 103 
s, 采用 本 文 所 建 IBSMC 方法 进行 20 次 双向 统计 过 
程 耗 时 增加 1.0 x 10-5 s, 仅 占 传统 蒙特 卡 罗 方 法 计 
算 耗 时 的 0.08%, 而 传统 蒙特 卡 罗 方 法 要 达到 相同 的 
倒 易 性 误差 精度 和 完整 性 精度 需要 的 模拟 光束 量 大 
约 为 每 个 单元 1 x 1019, 耗 时 为 3.9 x 104 s, 远 多 于 
IBSMC 方法 ， 

为 考核 所 建 IBSMC 方法 在 热 辐射 传输 问题 求 
解 上 的 优势 ， 选 取 如 图 4 所 示 的 模型 进行 计算 ， 模 
型 是 由 三 个 正方 形 表 面 组 成 的 腔 体 结构 ， 尺 寸 如 图 
所 示 ， 材 料 密度 为 500 kg/m3， 热 传导 系数 为 0.1 
W/(m-K)， 比 热 容 为 100 J/(kg:-K)。 外 部 边界 加 载 
500 KK 的 恒温 边界 条 件 。 辐 射 表面 采用 漫 灰 体 假 
设 , 表面 发 射 率 为 0.2, 每 个 表面 单元 模拟 光束 量 为 
1 x 103。 网 格 划 分 如 图 5 所 示 。 

采用 非 定常 求解 方法 求解 该 辐射 平衡 问题 ， 时 


间 方 向 采用 Newton-GMRES 隐 式 迭代 推进 求解 , 时 
间 步 长 取 0.01 s， 计 算 17000 个 时 间 步 。 计算 表 
明 ，16000 个 时 间 步 之 后 温度 场 不 再 变化 ， 说 明 问 
题 达 到 了 辐射 平衡 . 


1.0 
一 = 一 单元 1 
一 e- 一 单元 2 
0.5 一 一 单元 3 
了 一 +? 一早 元 4 
© 
x 
闵 0.0 幸 - 增 不 不- 不 让 一 增 不 一 夏 一 不 一 不 一 不 一 站 一 万 一 奢 一 矿 一 不 一 单一 
测 
MR 一 0.5 
一 1.0 
0 2 4 6 8 10 12 14 16 18 20 
双向 统计 迭代 次 数 
图 3 完整 性 误差 
Fig. 3 Error of the integrity 
9 
= 
> 


图 4 热 辐射 问题 计算 模型 
Fig. 4 Model for Radiative Heat Tansfer 


图 5 网 格 划 分 
Fig. 5 Mesh for Radiative Heat Tansfer 
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图 6 给 出 了 辐射 平衡 温度 分 布 随 双 向 统计 迭代 
步 数 的 变化 。(a)-(f) 分 别 代表 了 双向 统计 和 迭代 步 数 
为 0、1、5、10、15、 20 时 的 辐射 平衡 温度 分 布 . 理论 
辐射 平衡 温度 分 布 为 500 K 的 均匀 分 布 . 可 以 看 出 ， 
双向 统计 蒙特 卡 罗 方法 (n=1) 虽然 在 求解 精度 上 和 较 
传统 蒙特 卡 罗 方 法 有 明显 的 提高 ， 但 其 误差 还 是 相 
对 较 大 .而 和 迭代 双向 统计 蒙特 卡 罗 方 法 的 求解 精度 
随 迭 代步 数 的 增加 而 提高 ， 当 迭代 步 数 大 于 15 时 ， 
稳 态 温度 分 布 已 基本 接近 理论 值 . 


499 
—0.002 433:4 _0.002 
499.2 
498.8 _0.004 
> 498.6 
一 0.004 一 0.002 0 0.002 0.004 一 0.004 一 0.002 0 0.002 0.004 
X/m X/m 
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(e) (D 
图 6 不 同 迭 代 次 数 时 平衡 温度 场 分 布 


Fig. 6 Temperature distribution along the number of 


iterations 


平衡 温度 场 最 小 温度 与 理论 温度 之 差 、 最 大 温 

度 与 理论 温度 之 差 和 最 大 温度 与 最 小 温度 之 差 随 双 
向 统计 迭代 步 数 的 变化 如 图 7 所 示 , 由 图 7 可 以 看 
出 ， 当 双向 统计 和 迭代 步 数 为 1 时 ， 即 双向 统计 蒙特 
卡 罗 方 法 ， 温 度 差 虽 有 明显 下 降 ， 但 此 时 温度 差 仍 
然 相 对 较 高 。 随 着 双向 统计 和 迭代 步 数 的 增加 ， 温 度 
差 急 剧 下 降 ， 当 迭代 次 数 大 于 15 时 , 温度 差 趋 于 稳 
定 ， 已 基本 接近 于 零 ， 说 明 迭 代 双 向 统计 蒙特 卡 罗 
方法 可 以 有 效 提高 热 辐射 传输 问题 的 求解 精度 ， 且 
精度 较 双 向 统计 蒙特 卡 罗 方 法 有 较 大 提高 . 


; 一 一 最 小 温度 与 理论 温度 之 差 
1.5 一 一 最 大 温度 与 理论 温度 之 差 
一 一 最 大 温度 与 最 小 温度 之 差 
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0 2 4 6 8 10 
双向 统计 迭代 步 数 
图 7 不 同 迭 代 次 数 时 温度 误差 


Fig. 7 ‘Temperature distribution along the number of 


iterations 


3 结 论 


本 文 将 求解 热 辐射 传输 问题 的 双向 统计 蒙特 卡 
罗 方 法 分 解 为 归 一 化 过 程 和 双向 统计 过 程 ， 基 于 
迭代 思想 提出 了 和 友 代 双向 统计 楷 特 卡 罗 方 法 ， 分 
析 表 明 本 文 所 建 方法 可 以 进一步 提高 辐射 传递 因子 
的 倒 易 性 ， 从 而 以 较 小 的 耗 时 代价 进一步 提高 热 
辐射 计算 精度 . 通过 文中 算 例 和 分 析 可 以 得 到 以 下 
结论 

1) 本 文 所 建 迭代 双向 统计 蒙特 卡 罗 方 法 摆脱 了 
传统 双向 统计 蒙特 卡 罗 方 法 对 等 权 抽样 假设 的 依赖 ， 
具有 更 广 的 适用 性 。 

2) 迭代 双向 统计 蒙特 卡 罗 方 法 获得 的 辐射 传递 
因子 满足 完整 性 ， 其 倒 易 性 误差 的 其 级 随 双 向 统计 
迭代 次 数 的 增加 而 呈 线 性 下 降 趋 势 。 

3) 将 迭代 双向 统计 蒙特 卡 罗 方 法 应 用 到 含 空 
结构 的 辐射 /导热 厅 合 边界 问题 , 结果 表明 本 文 所 建 
方法 可 以 提高 热 辐 射 传输 问题 的 求解 精度 ， 并 且 求 
解 精度 随 着 双向 统计 和 迭代 次 数 的 增加 而 提高 . 
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